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SUMMARY 


A computational fluid dynamics (CFD) code with finite-rate 
reactions, FDNS2DR , has been developed to study the start 
transient of the Space Shuttle Main Engine (SSME) Fuel Preburner 
(FPB) operation. FDNS2DR is a versatile, pressure based, 
compressible/incompressible, implicit/explicit CFD code. It uses 
a time centered, time marching algorithm for time accuracy. An 
upwind scheme is employed for spatial discretization. The upwind 
scheme is based on second and fourth order central differencing 
with adaptive artificial dissipation. A state of the art two- 
equation turbulence model is employed for the turbulence 
calculation. A PARASOL chemistry algorithm is coupled with the 
point implicit procedure. In this study, FDNS2DR was benchmarked 
with three well documented experiments: a confined swirling 
coaxial jet, a non-reactive ramjet dump combustor and a reactive 
ramjet dump combustor. Excellent comparisons were obtained with 
the benchmark cases. The code was then used to study the start 
transient of an axisymmetric SSME fuel preburner. Excellent 
agreement was obtained for the temporal occurrence and magnitude 
of the temperature spikes in comparison with results obtained 
from an Aerojet experiment. Furthermore, it was also found in 
this study that an appreciable amount of unburned oxygen entered 
the turbine stages. 

The significance of this study is that a time accurate CFD 
design tool has been developed to accurately predict the severe 
thermal gradients which are impressed upon the fuel preburner and 
turbine during the start transient. Previously, various factors 
were suspected to cause the ignition temperature spikes, such as 
the collection of unreacted propellants, a mismatch of the fuel 
and oxidizer supply rate, and mixture maldistribution prior and 
during ignition. These factors can now be studied, 
quantitatively. Numerical experiments can be designed for study 
with FDNS2DR to alter the ignition and shut down characteristics 
to attenuate the possible thermal loads. 
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NOMENCLATURE 


A a NxN species jacobian matrix 
A^ chemical symbol of the 1th species 
A r Arrhenius constant 

ASI subscript denotes augumented spark igniter flow properties 
B a N element column vector 

C-n heat capacity for species n 

Cj turbulence model constant, *1.15 for k-e{T> model 

*1.43 for standard k-e model 
C 2 turbulence model constant, *1.90 for k-e{T} model 

*1.92 for standard k-e model 
C 3 turbulence model constant, *0.25 

C 4 turbulence model constant, =0.6 

turbulence model constant, *0.09 
c local speed of sound 

D outer pipe diameter of ramjet dump combustor 
d dissipation term 

E activation energy 

F convection and diffusion fluxes 

G geometrical matrices 

h enthalpy 

INJ subscript denotes faceplate injector properties 
J Jacobian of coordinate transformation 
J n diffusion flux for species n 

K chemical reaction rate 

k turbulent kinetic energy 

M molecular weight 

N total number of species 

P Pade' approximants in ^th order 

p static pressure 

p c fuel preburner chamber static pressure 
Pr turbulent kinetic energy production 
Q Pade' approximants in 5th order 
R gas constant 

REC adaptive dissipation parameter 

S fuel preburner inlet area 

Sg source term for equation q 

t static temperature 

T* ratio of static to room temperature 
t time 

U contravariant velocity 

u x-direction mean velocity 

v y-direction mean velocity 

u' updated transient inlet x-direction mean velocity 

v' updated transient inlet y-direction mean velocity 

W n mass production rate for species n 

W A mass flow rate of either fuel or oxidizer 

W T total mass flow rate 

x physical coordinate 

Y a/M 

G energy dissipation function 


v 


-< ■a> ■© -o 


(,7j computational coordinate 
V e ff effective viscosity 
V- molecular viscosity 

turbulent eddy viscosity 
e turbulent kinetic energy dissipation rate 

°q turbulence model constant for equation q. (a u = 1.0, o -*0.9, 
<7^=1. 0, <7^=0.8927 and ae=1.15 for k-e{T} model, 

<7^=1. 0 and a £ = 1.3 for standard k-e model 
o species mass fraction 
density 

updated transient inlet density 
pressure relaxation parameter, =10.0 
temperature dependency for Arrhenius equation 
r ( chemistry time integration step 
^it, stoichiometric coefficients 

{ } denote functionality 
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INTRODUCTION 


Numerous flow anomalies have been found in both the fuel and 
oxidizer preburners and turbine hot end parts. Cracks have been 
found on housings, Kaiser-hat nuts, struts, sheetmetal, nozzles 
and shrouds. After many studies and modifications to the 
preburner components, several problems remain with the SSME 
preburner operation, such as cracking of the first and second 
stage blade shanks, of the blade root fillets and under the blade 
platform. The deunage is more severe on the fuel side since it 
operates at a higher power level and temperature. It is 
postulated that the extreme thermal environment during the start 
up transient produces the initial blade cracks; the severity of 
the cracks will depend on the specific transient operations the 
engine experiences. The cracks can continue to grow due to both 
high cycle fatigue and low cycle fatigue as the number of turbine 
load cycles accumulate. Other flow anomaly theories such as 
inefficient atomization and possible unmixedness of the gaseous 
oxygen during steady state operation are considered to be of 
secondary importance in creating locally excessive heating 
problems . 

Computational analyses of the axisymmetric and three- 
dimensional FPB flowfield (ref. 1) indicated that under steady 
operating conditions, striated flows go through the preburner 
chamber with little turbulent mixing between the parallel shear 
layers. Mixing at the injector face is thought to be very rapid, 
and assuming that some mixing occurs immediately allows the 
entire flowfield to be simulated. These observations and 
assumptions are consistent with a computational analysis of a 
single injector element (ref. 2), in which the percentage of 
oxidizer atomization, vaporization and fuel/oxidizer reaction 
reached 100 percent almost immediately below the injection inlet. 
Results from such analyses show that striations severe enough to 
cause local streamlines to be at near stoichiometric mixture 
ratios should not be present at the inlet to the turbopump under 
steady-state operation of the preburner. Temperature 
measurements at the turbine inlet (ref. 3) shown in Fig. 1 show 
that two temperature spikes exist under start-up. Estimates of 
thermocouple response suggest that the actual temperatures are 
sufficiently high to create blade cracking (ref. 4). A 
transient fuel preburner thermal flowfield computational analysis 
was performed in this study in order to understand and make 
recommendations for improving the component life of the SSME 
engine. 
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TRANSIENT FUEL PREBURNER OPERATION 


The configuration of the FPB is shown in Fig. 2. The flow 
passage from the fuel preburner oxidizer valve (FPOV) into the 
oxidizer manifold above the injector faceplate is shown in this 
figure. The oxidizer flow to the augmented spark igniter (ASI) 
is supplied from a bleed line from the FPOV. The FPOV is shown 
in Fig. 3. The valves which control the engine flows are shown 
schematically in Fig. 4. Notice there is not a fuel control 
valve for each preburner. Flows from the fuel and oxidizer 
manifolds enter the FPB chamber through the injector faceplate, 
the layout of which is shown in Fig. 5. These figures are from 
references 3 and 5. 

The coaxial injectors and the ASI are designed to be high 
intensity mixers. The fact that additional hydrogen jets were 
required to mix the ASI flow indicate that the original design 
was not as sound as anticipated, nevertheless the current 
configuration does apparently avoid excessive hot spots on the 
turbopump dome. When simulating the preburner flowfield, it is 
not computationally feasible to resolve each flow injector or the 
detailed convective flows within the ASI. Rocketdyne has 
reported (ref. 6) experimental measurements of temperature, 
velocity, density gradients, and OH concentrations for single 
coaxial elements and for groups of four coaxial injector flows. 
These measurements showed that injectors create back-flow 
recirculation up to the injector faceplate, with the hydrogen 
from the bleed holes constituting the bulk of the fluid which is 
recirculated. If unmixed lumps of oxygen were present, the hot 
spheres of radiating OH shall be visible as the hot balls of 
gases react. Such flow structures were not observed, hence, 
large unmixedness effects of the oxygen probably does not exist. 
Ignition probes provide sparks that cause ignition. Initial 
flames move around until recirculating flows cause the entire 
near face flow to be ignited. The Injector faceplate then acts 
as a flameholder. Incomplete mixing along the centerline of the 
coaxial jets cause the flow to be cool until sufficient hydrogen 
diffuses in to change the mixture ratios. These cool regions are 
roughly located from two to six Inches downstream of the 
faceplate. These experiments support the contention that 
unmixedness is not a major factor in preburner flows. Detailed 
experiments of ASI flows are not available, but the ASI was 
designed so that tangential velocities of hydrogen cause the flow 
to be intensely mixed. Assuming that the hot spots on the 
turbopump dome were caused by ASI flow, the extra hydrogen which 
would have to remain unmixed suggests that a hot core surrounded 
by very cold, almost pure, hydrogen would have to exist. The 
modification which consisted of adding extra hydrogen jets on the 
baffle tips was apparently adequate to solve this problem, 
although, the simultaneous change in ASI start-up procedure may 
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Fig. 5 The layout of SSME fuel preburner injector faceplate. 
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have been the actual cause of reducing the thermal loads In the 
turbopump dome. Computational analysis for single coaxial 
elements with the ARICC code (ref. 2) also suggests that 
atomization and combustion is complete in about two inches, 
although radial concentration gradients may still exist out to 
several more inches. An alternative analysis may also be 
considered by using Rosner's treatment (ref. 7) of supercritical 
combustion. This more elementary model treats the process as a 
one-dimensional mass diffusion with laminar diffusion 
coefficients. For a mean drop (or lump) size of 150 ftm, 
combustion should be complete in about two inches. These 
analyses are admittedly approximate but they offer no 
justification for arbitrarily delaying the mixing (by fixing 
transport coefficients) until oxygen globules are present at the 
preburner exit, during steady state operation. Since unmixedness 
appears to unlikely cause of turbine blade heating, what other 
explanations can be postulated? This is the motivation for a 
transient FPB study. 

The turbine inlet temperature measurements reported in ref. 

3 and the one-dimensional, transient start-up analyses performed 
by ref. 8 indicate that transient flow phenomena may be the cause 
of high thermal loads on turbine blades. Operations of the 
preburners are designed to have a fuel lead, but the starting 
flows of oxygen may be much different from those at steady state 
operation. Two scenarios can be used to explain the transient 
flow phenomena of the fuel preburner operation. The first 
scenario recognizes that under both steady and transient start-up 
the hydrogen is at gaseous state. But gaseous oxygen (GOX) is 
probably the first oxidizer to flow into the preburner although 
LOX is injected at steady-state. Actually, these flows are even 
more complicated because valve sequencing and transient heat 
transfer must be considered in detail to describe these multi- 
dimensional, unsteady flows. Interpreting these flows is even 
further complicated by the fluid state being supercritical with 
respect to pressure. If LOX flashes to GOX as it flows into the 
unchilled LOX dome, then GOX will flow through the injectors for 
some brief period of time. Such GOX flow might make the 
preburner leaner than expected before steady state flow is 
established. The second scenario considers the response of the 
flow in the fuel supply line to the increase of pressure in the 
combustion chamber. Remember there is no fuel control valve in 
the fuel supply line, as shown in Fig. 4. The engine main fuel 
valve begins to ramp open at the start command, and the FPB 
oxidizer valve begins to open at about 0.1 of a second. 
Propellants begin to accumulate in the preburner but do not 
react. Combustion does not actually begin until about 0.4 of a 
second, when the turbine inlet temperature begins to rise. The 
preburner chamber pressure also starts to rise after about 0.4 
seconds during the priming process. The pressure rise rate is 
quite steep after the priming event. A pressure wave would be 
felt upstream of the faceplate. The fuel supply line would feel 
the pressure increase and shut down the fuel supply momentarily. 
To compensate for the fuel flow loss, the oxidizer flow rate will 
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have a temporary increase. This combination of fuel flow loss and 
oxidizer flow gain would create a temporarily leaner mixture 
ratio in the chamber. 

Aerojet (ref. 3) studied the engine test data and concluded 
that the two temperature spikes which are observed on start-up 
are due to burning of propellants accumulated in the combustion 
chamber and extinguishment caused by cold slugs of fuel flowing 
into the FPB before the steady-state is achieved. Seymour (ref. 
8) has made an extensive DTM study and concluded that temperature 
spikes are closely associated with the decline of the fuel supply 
which occurs during start-up. These observations are thought to 
favor the second scenario, although they do not totally exclude 
the possibility of the first one. The DTM simulation also 
provides the precise sequencing of control valves and the 
reaction of the fuel supply line to the pressure wave, which are 
thought to be the dominating factor of the transient FPB 
operation by the second scenario. It is therefore proposed in 
this study that the result of the DTM simulation be used as the 
input for the transient FPB thermal flowfield study. 

Since the results of Seymour's study were not available at 
the beginning of this investigation, an analysis using the first 
scenario as the possible explanation of the transient mechanism 
was initiated. Two codes, COBRA/TRAC and FDNS were evaluated for 
use in describing the transient flow of GOX into the inlet 
manifold to the preburner and the subsequent chilling of the 
manifold. It was determined that either of these codes could be 
used for making the required analysis, if suitable oxygen 
property data were available. A general procedure for the 
calculation of the thermodynamic and transport properties for two 
phase oxygen was prepared and is included in the Appendix for 
future usage. However, when Seymour's analysis revealed that the 
relative flows of the two propellants could possibly explain the 
temperature spikes, development of the first scenario was 
terminated in order to evaluate the second one. 

Based on the considerations just presented, it was 
determined that the transient fuel preburner could best be 
simulated by including finite-rate combustion reactions in the 
FDNS code and using the DTM to represent the unsteady upstream 
boundary conditions. The reacting FDNS code was developed, 
validated with several test cases, and the transient operation of 
the fuel preburner was simulated. The results of these analyses 
are described in the following sections of this report. 
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FDNS METHODOLOGY 


Governing Equations 

The basic equations employed in this study to describe a 
general combustor flowfield are the axlsymmetrlc, multi -component 
conservation equations. A generalized form of these equations 
written in curvilinear coordinates is given by 

( 1/J) ( dpq/dt) - 3[-pU iq + PeffGijOq/aejn/afi + S q (1) 

where J, and G^j represent the Jacobian of the coordinate 
transformation, contravariant velocities and diffusion metrics 
respectively. They are written as: 

J = d((,v)/d(x,y) 

U A = (uj/J) (a^/axj) 

G ij * Oi 1 /ax k ) (a^j/ax k )/j ( 2 ) 

f*eff = + P t^ a q is the effective viscosity when the 

turbulent eddy viscosity concept is employed to model the 
turbulent flows. « pC^k 2 /e, is the turbulence eddy viscosity 

and and a„ denote turbulence modeling constants. Source terms 
Sg are given^by 


S q = (1/J) 


0 

-P x + ^eff( u j>xJ 
-Py + *(P e ff 

Dp/Dt + * + ZJ n C pn vT - Ih n W n 

P(P r “ *) 

p(e/k) (C 1 P r -C 2 e+C 3 T* C4 Pr 2 /e) 
w n 

n = 1 , . . . . ,N 


( 3 ) 


The equation of state for an ideal gas is used to close the 
above system of equations. 


Numerical Schemes 


The spatial accuracy of numerical schemes for finite 
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difference or finite volume Navler-Stokes solvers applied to 
convection dominated flow problems depends strongly on the order 
of approximation in discretizing the convection fluxes. In 
recent years, many numerical schemes have been developed for high 
to low speed flows based on the concept of upwind differencing. 
The flux vector splitting schemes (refs. 9 and 10), the flux 
difference schemes (refs. 11 - 13) and central plus artificial 
dissipation schemes (refs. 14 - 16) are the most commonly used. 
All these schemes approximately follow the characteristics of the 
flow by an upwinding treatment of the numerical fluxes, while 
attempting to satisfy the conservation laws. Numerical accuracy 
of second-order or higher can be attained throughout the entire 
flow field except near regions of sharp gradients, e.g. shock 
wave discontinuities. It is also a common practice to use first- 
order upwind schemes near shock waves to maintain monotonicity of 
the solutions; this practice is justified by the total variation 
diminishing concept (ref. 13). 

It is also known that both flux splitting and flux 
difference methods require certain matrix operations in order 
to construct an upwind scheme based on the eigenvalues of the 
flux Jacobian matrices. On the other hand, the methods of 
MacCormack and Jameson require no extra matrix operations. 
Additional savings in computational time are achieved when 
artificial dissipation methods are employed. As has been pointed 
out by Pulliam (ref. 17), artificial dissipation methods can be 
very accurate and do not exhibit numerical oscillations near 
shock waves, if they are constructed to reflect the same features 
as the flux vector upwind methods. 

In the present development, an upwind scheme was employed to 
approximate the convective terms of the momentum, energy and 
continuity equations; the scheme is based on second and fourth 
order central differencing with artificial dissipation. First 
order upwinding is used for the species and turbulence equations, 
since the parameters involved are positive quantities. Different 
eigenvalues are used for weighing the dissipation terms depending 
on the conserved quantity being evaluated, in order to give 
correct diffusion fluxes near wall boundaries. This procedure is 
different from those proposed in other works (refs. 13 - 15) in 
which the sum of the absolute value of the convection velocity 
and the local speed of sound is used to weigh the dissipation 
terms. For simplicity, let us consider fluxes in (-direction 
only. That is, 

3F/9( s (F^ + j- ~ (&i+i/2~ ^i-1/2^ (*) 

The dissipation terms are constructed such that a fourth- 
order central and fourth-order damping scheme is activated in 
smooth regions, and a second-order central and second-order 
damping scheme is used near shock waves. Since the Jacobian 
matrices of the Euler fluxes have eigenvalues of U, U+c and U-c, 
it is sufficient to use the magnitudes of these eigenvalues to 
weigh the dissipation terms to maintain the smoothness of the 
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solution without losing accuracy. It is suggested in this study 
that |U|+c should be used for the continuity equation and ]U| 
should be used for other transport equations. A general form of 
the dissipation term is given below. 

d i+l/2 “ 0 . 25 £ 2e j | pU | + e 2 ( |U|+c) ] 1+1/2 
(Si+i " ^i) + [fiad-eiJMAXto.s 
Asp ( | U | + | V | ) , 2 | pU | ] + £ 4 ( |U|+C) ] i+1/2 


(qi -1 - 3q A + 3q 1+1 - q 1+2 ) (5) 

where a ± > |p 1+1 - 2p A + Pi_ 1 |/(P i+1 + 2p A + Pi_ 2 ) 
y i+l/2 “ MAX(a2» a^ + j) 
d 2 = MAX[REC,MIN{1 .0,25y i+lr / 2 ) ] 
d 4 = MAX (0.0, 0.01-0.251)2 +1/2 ) 

Different values for e 2 , e 3 and e 4 are used for the 

continuity, energy and momentum transport equations. Table 1 
summarizes these parameters. 


Table 1. Dissipation Parameters 



Momentum & Energy 

Continuity 

£ 1 

d l 

0 

£ 2 

0 

u i+l/2 

e 3 

0.015 

0 

e 4 

0 

d 4 


A pressure based solution method was selected so that a wide 
range of flow speeds could be analyzed with the same code. 
Successful results of viscous flow computations using pressure 
based methods have been reported (refs. 18 and 19) • For high 
speed flow cases, a hyperbolic pressure correction equation was 
employed by perturbing the density in the mass conservation 
equation. This provides a smooth transition from low to high 
speed flow characteristics. For time accuracy, a time-centered, 
time-marching scheme with a multiple pressure corrector algorithm 
was employed. In general, a noniterative time-marching scheme 
was used for time dependent flow computations; however, 
subiterations can be used if necessary. 

To solve the system of nonlinear coupled partial 
differential equations, Eq. 1, finite difference approximations 
were used to establish a system of linearized algebraic 
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equations. A relaxation solution procedure was employed to 
couple the governing equations. For convenience, Eq. 1 was 
rewritten as: 


(1/J) (3pq/3t) - -3F 1 /8^ 1 + S q * R q (6) 

First, Eq. (6) is discretized in time with a time-centered 
(Crank-Nicholson) scheme. That is, 

(l/JAt)[(pq) n+1 - ( pq) n ] = (R q n+1 + R q n )/2 (7) 

where the superscript n denotes the current time level. If a 
sub-iteration procedure within a time step is applied, the 
following linearization can be incorporated. 

(pq) n+1 = (pq) k + p n Aq k (8) 

R q n+1 = (3R q /3q) k Aq + R q k (9) 


where the superscript ^ denotes the k-th sub-iteration. With the 
above approximations, the final form of the time-marching scheme 
can be written as: 

[(p n /JAt) - ( 3R q /3q) k ] Aq k 

* — (1/JAt) [ (pq) k — ( pq) n ] + (R q k + R q n )/2 (10) 

The solution at time level n+1 is then updated by: 

q n+l = q k+l = q k + Aq k (11) 

When k = 1 is selected, a non-iterative time-marching scheme is 
used. As reported in ref. 18, the non-iterative time-marching 
scheme with a multi-corrector solution method can provide time 
accurate solutions for transient flow problems. This multi- 
corrector procedure will be described below. 

A simplified momentum equation was combined with the 
continuity equation to form a pressure correction equation. This 
pressure correction equation exhibits elliptic behavior for low 
speed flow and becomes continuously more hyperbolic as flow speed 
Increases. The simplified momentum equation can be written as: 

3pu i /3t = - vp * (12) 

or, in discrete form, 

u'i a -j3(At/p)vp' (13) 

where 0 represent a pressure relaxation parameter (p * 10 was 
used for all test cases in this paper) . The velocity and density 
fields in the continuity equation are then perturbed to form a 
correction equation. That is, 


13 



V(pu i ) n+1 = 7[{p n + P')(u n i + u’i)] = 0 (14) 

By neglecting the p'u' terms, the following equation results. 

*(u ip') + vjpu'i) = - v(p Ui ) n (15) 

Substituting Eq. (13) into Eq. (15) and letting p*« p'/RT, the 
following pressure correction equation is obtained. 

(u i /RT)p' ] - v(j9At vp') « - v(p Ui ) n (16) 

To provide smooth shock solutions the adaptive dissipation terms 
described above were added to the right hand side of Eq. (16). 
Once solution of Eq. (16) is obtained, the velocity field and the 
pressure field are updated through Eq. (13) and the following 
relation. 

p n+l = p n + p , (17) 

The density field is then updated by applying the equation of 
state. To ensure that the updated velocity, density and pressure 
fields satisfy the continuity equation, the above pressure 
correction solution procedure is repeated several times (usually 
4 times are sufficient) before marching to the next time step. 
This represents a multi-corrector solution procedure. 

To solve the system of linear algebraic equations, a semi- 
implicit matrix solver (ref. 20) was utilized. This matrix 
solution algorithm is a modified Stone's method (ref. 21). This 
method has been reported to be very efficient for matrix systems 
resulting from complex curvilinear coordinates systems including 
high aspect ratios and high grid skewness. However, this method 
does require a large amount of memory. 


PARASOL, a Finite Rate Chemistry Algorithm 

The finite rate chemistry source terms were evaluated with a 
point implicit procedure before the species equations were 
solved. The source term was solved by an improved algorithm, 
PARASOL (Pade' Rational Solution). This algorithm was first 
proposed by ref. 22 for stiff chemistry systems. The improvement 
made in this study was to incorporate LINPACK (ref. 23) , instead 
of the original Crout's method, to solve the matrix system. A 
benchmark case involving twelve reactions and eleven species in a 
streamtube calculation resulted in fifty times reduction in CPU 
time in comparison to the Crout's method. 

Consider a general chemical system of N species undergoing m 
simultaneous elementary reactions 
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k=l,2. 


m 


(18) 


N 

z 

1=1 


y ik A i 


N 

z 

1=1 


u ik k i 


The source term S„4 for volumetric chemical reaction mass 
production rate of species i is 

Sgi = pfdai/dt) 

m N 

= M a l (I'ifc-i'ik) [K fk D (py ) y rk 
k=l r=l 


N 

- K bk n (py r ) y rk ], (19) 

r=l 

1=1,2, N 

Eg. (19) is a set of coupled, nonlinear system of stiff 
ordinary differential equations for general combustion systems. 
The stiffness comes from the inherent, widely disparate time 
constants of the nonlinear terms in Eq. (19). The stiffness can 
also be interpreted as the drastic change of local eigenvalues 
occurred near a flame front (ref. 24). PARASOL started by casting 
Eq. (19) in a general form 

dy/dt = f{y, T, p) (20) 

Nonlinear instability will occur if Eq. (20) were to be 
solved directly with even modest time steps. Locally linearizing 
Eq. (20) through the first order terms, and asstiming the 
temperature change is small within an integration step, yields a 
system of linear ordinary differential equations with constant 
coefficients and can be regrouped as 

dy/dt = A y{ t } + B (21) 

The formal solution of Eq.(21) at t=r is 

y{r} = exp{ rA) (y Q + A _1 B) - A -1 B (22) 

In Eq . ( 22 ) , the matrix exponential term can not be evaluated 
directly but may be approximated as polynomials, and the form of 
the approximation will determine whether it is accurate, 
efficient and stable. For a family of unconditionally stable 
schemes, the Pade' approxlmants are found to be more efficient 
than the other approxlmants. The Pade' approximation of the 
matrix exponential of Eq.(23) can be expressed as 

exp(rA) = Q _1 P + 0(rA)^ +a+1 (23) 

while 
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X (X + <5 - k) ! 

P - l ( rA) k (24) 

k=0 ( X + 6) !k! (<5 - k) ! 

and 

S (X + 6 - k)'.6l 

Q = I (-rA) k (25) 

k=0 (X + a) !k ! (6 - k) ! 


The truncation error means that these formula must agree with the 
exponential power series for at least X + 6 + 1 terms. 
Substituting the approximants into Eq. (22), the integration 
formula becomes 


y(h) = Q _1 [Py 0 +MP - Q ) A -1 B ] 


(26) 


The diagonal Pade ' approximants (6 - X) have the smallest 
truncation error of these approximants and are A-stable. The 
integration formula Eq. (26) can be derived for any order of 
truncation. Fifth order polynomials (p • q = 2) were found to be 
sufficient in this study. An automatic species Jacobian writer 
was also developed in this study. 
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Turbulence Modeling 


k-e(T) Turbulence Model 

The two-equation k-e turbulence modeling approach was 
employed In the present study to describe the flowfield. The 
conventional wall function approach was used for near wall 
boundary conditions. It has been shown that in the extended k-£ 
model (ref. 25), two time scales (the production range time 
scale, k/Pr, and the dissipation rate time scale, k/e) can be 
Included in the dissipation rate equation to cause the 
dissipation rate to respond more effectively to the mean strain 
than it does in the standard k-e model. The net effect of that 
energy transfer function is to enhance the development of e when 
the mean strain is strong and suppress its production when the 
main strain is weak. The source term for the dissipation rate 
equation is: 

S £ = p(C 1 Pre/k - C 2 e 2 /k + C 3 Pr 2 /k) (27) 

The last term represents the energy transfer from large 
scale turbulence to small scale turbulence controlled by the 
production range scale and the dissipation rate time scale. This 
formulation enables the dissipation rate to control the 
development of the turbulent kinetic energy more effectively. It 
has proved to be effective for isothermal flows such as boundary 
layer type flows and general elliptic type flows. For non- 
isothermal flows such as a chemically reacting flow in a ramjet 
dump combustor configuration, the isothermal standard k-e and 
extended k-e models predicted too rapid a decay of the centerline 
kinetic energy profile. The extended k-e model did predict a 
slower centerline velocity decay than the standard k-e model. It 
was therefore postulated that a temperature dependent term could 
be added to the last term of Eq. (27) for non-isothermal flows, 
without changing the modeling constants already determined for 
the isothermal turbulent flows. The value of the temperature 
dependency term would be unity for isothermal flows and provides 
no effect for the turbulent flowfield. The rationale behind this 
new turbulence model (k-e(T)) is that the temperature rise in the 
reacting flow would cause an increase in the dissipation rate, 
and the net effect would be a reduced effective viscosity and a 
longer recirculation zone in a recirculating flow environment. 

The final expression for the dissipation rate equation source 
term is 

S e = p ( CjPre/k - C 2 e 2 /k + C 3 (T*) C4 Pr 2 /k) (28) 

The final values of these model constants for k-e(T) 
turbulence model are: 02=1.15, C 2 =1.9, C 3 =0.25 and C*=0.6. 

Note that the temperature dependence has been established only 
for the reactive ramjet dump combustor experiment. 
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Chemistry Modeling 


One Step Reversible H2/02 Kinetics 

To study the reactive hydrogen/oxygen system with the 
elliptic Navier-Stokes equations, the introduction of a detailed 
chemical reaction mechanism is often cost prohibitive. 
Simplification of the reaction system while maintaining Important 
combustion characteristics is a rational strategy for reactive 
flow modeling. This idea has been developed as a quasi -global 
kinetics model for complex hydrocarbon fuels like NO. 2 fuel oil 
and coal derived fuels (refs. 26 and 27). For hydrogen/air 
combustion, a seven reaction mechanism has been proposed in ref. 
28. Recently a two reaction global (ref. 29) was proposed for 
supersonic Scramjet combustor calculation by comparing the 
ignition delay times with a 28 reaction mechanism. This two 
reaction model consists of four species, therefore five species 
equations are needed for a hydrogen/air combustion calculation. 

In this study, a one step reversible H2/02 kinetics model was 
proposed and the reaction model is given by 

H 2 + 0.5 0 2 = H 2 0 (29) 

mhe reaction rates are functions of temperature and are 
usually expressed in a general Arrhenius form. For example, the 
forward reaction rate can be expressed as 

K f = A r T^exp{-E/RT) (30) 

The rate constants for the forward reactions are: 

A r =l . 1068E6 , r=0.5 and E=8568.9; and for the reverse reaction: 
A r =1.0E5 and E=19870.0. 

The rate constants were obtained through parameter 
optimization by comparing the ignition delay times and the 
temperature histories with those predicted from a 28 reaction 
mechanism. Fig. 6 and 7 shows the comparison of the one-step 
model with the 28 reaction mechanism for H 2 -air combustion. The 
predicted temperature profiles and the Ignition delay times for 
the one step model agreed well with those from the 28 reaction 
mechanism. It is not implied that these simplified kinetics rate 
constants have general applicability, merely that they are useful 
under conditions where they have been verified. 
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Fig. 6 Comparison of temperature histories for 
H2~air combustion at 1 atm. 
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Fig. 7 Comparison of temperature histories for 
H2~air combustion at 1 atm. 
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CODE VERIFICATION 


A Confined Swirling Coaxial Jet 


The ASI unit of the SSME FPB is a confined, small injector 
and combustion chamber. It initiates ignition of propellants in 
the suddenly expanded FPB combustion chamber. The injector has a 
single pair of impinging oxidizer orifices and eight hydrogen 
orifices directed tangentially around the oxidizer. The 
tangential motion of the hydrogen flow makes the Igniter flow a 
confined suddenly expanded swirling jet. Roback and Johnson 
(ref. 30) have accumulated extensive experimental data for 
swirling turbulent flows in confined suddenly expanded coaxial 
jets. Their experiment is therefore chosen for validation of 
FDNS2DR with respect to a confined swirling jet. 


Experiment 

Roback and Johnson (ref. 30) documented the flowfields for 
the configuration as shown in Fig. 8. Experimentally, the 
swirling flow channel consisted of two coaxial inlet pipes with 
the swirling guide vanes installed between the inner and outer 
pipes. The inner and outer pipes had radius of 12.5 mm and 29.5 
mm respectively. The inlet channel was followed by a sudden pipe 
expansion with expansion ratio at about 1:2. The radius of the 
downstream pipe, R 0 , was 61mm. The flow chosen for documentation 
has Reynolds numbers of 15,900 and 47,500 for the inner and 
annular streams respectively. These values are within the 
turbulent flow range and are typical of gas turbine combustors. 
For these experiments, water was the working fluid. The use of 
water limits the simulation in that it eliminates combustion and 
multiple species as variables to be considered. In this study, 
this experiment was only used to justify the ability of FDNS2DR 
to compute swirl. Multiple species mixing and combustion were 
evaluated with the dump combustor experiment. Flow condition 1 
of the documentation was chosen for the code evaluation. 


Results and Discussion 

In order to avoid geometrical and flow complexities upstream 
of the expansion plane, where a three-dimensional flow field was 
expected and detailed experimental data were not available, the 
present computation employed the first data plane, which is 5 mm 
downstream of the expansion, as the inlet boundary where detailed 
experimental data were provided (ref. 30). A grid size of 61x41 
was used for this test case. The grid layout is shown in Fig. 9. 
The exit boundary was located 1800 mm downs tream of the 
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expansion. Notice the ratio of length to height of the grid 
layout in Fig. 9 was not plotted to the actual scale. The grid 
points were clustered near the inlet boundary where rapid flow 
development was expected. 

The inlet swirling velocity generated by the swirling guide 
vane created a central recirculation zone along the pipe center 
line downstream of the expansion plane. This central 
recirculation zone was accompanied by a corner recirculation 
region downstream of the step. Comparisons of the development of 
the center line velocity are given in Fig. 10. The extended k-e 
model predicted much better results than that of the standard k-e 
model. Figs. 11 and 12 show the comparison of mean axial 
velocity and the mean azimuthal velocity at various axial 
stations, respectively. The extended k-e model predictions again 
performed better than those of the standard k-e model. The 
ability to predict the low wave number turbulence is believed to 
be the reason why the extended k-e model outperform the standard 
k-e model. 




Fig. 8 Sketches of confined swirling coaxial jets test section 
Inlet region with velocity and coordinate system. 
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Fig. 10 Comparison of mean axial Telocity along test section 
centerline, symbols, experiment; solid line, k~£{T( model; 
dotted line, standard k~£ model. 


25 



Z - 102MM 




M/S 





A Non-reactlve Hydrogen/Alr Ramjet Dump Combustor 


The ability to accurately predict the multiple species 
mixing and finite rate combustion of FDNS2DR are demonstrated in 
two test cases: a non-reactive hydrogen/air ramjet dump 
combustor, and a reactive hydrogen/air ramjet dump combustor. A 
tremendous effort has been devoted to the development of Improved 
numerical techniques for design and analyzing such systems in 
recent years. This is mainly due to the Increased demand for 
higher performance and the tighter system size and weight 
constraints on future propulsion and combustion systems in 
general, and of the National Aerospace Plane in particular. A 
critical requirement in achieving the future demand is the 
ability to accurately describe the combustor flowfleld for varied 
design parameters. Unfortunately, the physical phenomena 
occurring in a ramjet dump combustor are quite complex and 
present quite a task for numerical tools and computer resources. 
It is for this reason that modular models (refs. 31 and 32) were 
first used to predict ramjet combustor flowfields. Such models 
assumed that the entire flowfleld can be broken down into several 
neat elements: the recirculation zone, the core region, the shear 
layer and the mixing region. The recirculation zone can then be 
represented by a perfectly stirred reactor to stabilize the 
flame, and the central core can be described by parabolic 
equations. These two regions can then be joined by a shear layer 
model and the parabolic mixing model to complete the calculation. 
The modular model has obviously utilized assumptions regarding 
the structure of the flowfleld, e.g., the size of the 
recirculation region. It will be shown later in this study, that 
the recirculation length might vary by a factor of two, depending 
on whether the flow is reactive or non-reactive. Nevertheless, 
the modular approach is a computational simple and physically 
sound model which is useful in initial parametric analyses. 

More complete models that recognize the spatially elliptic 
nature of ramjet dump combustor flowfields require the solution 
of the full Navier-Stokes equations. These models describe the 
entire dump combustor flowfield without resorting to the modular 
components. An excellent pioneering work was reported by 
Drummond (ref. 33), his analysis of the dump combustor flowfield 
Included the solution of the spatially elliptic and temporally 
parabolic form of the governing equations. An unsplit MacCormack 
finite difference scheme (ref. 34} was used to Integrate the 
system equations with time stepping. An algebraic eddy-viscosity 
model was used for the turbulence, and a decoupled, instantaneous 
reaction was utilized for the chemical reaction. The results 
from that study predicted the flowfield to be somewhat more 
diffusive than observed in the experiments. Drummond concluded 
that the overprediction of the mixing rate and large diffusion 
were caused by the turbulence model. Trinh (ref. 35) has 
performed a similar analysis on a non-reactive dump combustor 
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simulation. A standard two-equation k-e turbulence model (ref. 
18) and a staggered grid finite difference solver (ref. 19) were 
used to perform the calculation. Improved mixing rates and 
recirculation region size were obtained with the more 
sophisticated k-e turbulence model, but the reacting cases were 
not analyzed. 


Experiment 

A reunjet dump combustor simulation experiment was carried 
out at the Arnold Engineering Development Center (ref. 36). The 
dump combustor simulation facility is shown in Fig. 13. Air was 
injected from the inner 2.07 in. diameter pipe at a mass flow 
rate of 0.58 lbm/sec. The inner air nozzle was surrounded by a 
5.24 in. diameter outer pipe from which gaseous hydrogen was 
injected at a mass flow rate of 0.002 lbm/sec. The one- 
dimensional velocity of the primary and secondary stream was 335 
ft/sec and 3 ft/sec, respectively and the inlet temperatures were 
nominally 525 deg. R. for the nonreactive case and 1300 deg. R. 
for the reactive case. The nominal static pressure in the 
exhaust duct was 13.7 psia. Measurements were made at axial 
stations from zero to six duct diameters from the primary nozzle 
exit plane. The diameter ratio of this configuration was 2.5. 


Results of the Non-reactive Dump Combustor Calculation 

The analysis of the non-reacting dump combustor flowfield 
was carried out in a physical domain, beginning one half outer 
pipe diameter upstream of the step and concluding six diameters 
downstream of the step. The central primary flow was air, and 
the secondary flow coming out of the step was hydrogen. Fig. 14 
shows the final computational grid (60x31) selected from a series 
of grid dependency tests. The grid lines were clustered near the 
wall and the corner of the step to resolve the boundary layers. 
The predicted non-reacting dump combustor steady state vector 
field is shown in Fig. 15. Notice the recirculation region 
downstream of the secondary flow step. The convergence history 
is given in Fig. 16. Approximate convergence was obtained after 
the original value of the residual had dropped four orders of 
magnitude. It took 525 CPU seconds and 1066 Iterations for a 
typical k-e(T) turbulence model calculation on a CRAY XMP 
computer. The secondary flow velocities were estimated for this 
simulation, since negative velocities were reported (ref. 36). 
Estimates were based on the overall secondary mass flow rates. 
This was also the case for the reactive flow calculation. 

Fig. 17 shows the comparison of the experimental and the 
computed centerline velocities. The temperature dependency of 
the k-e(T) model was negligible on the non-reacting case since 
the flow was at room temperature. It therefore behaves like the 
extended k-e turbulence model. Both turbulence models predicted 
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Fig. 13 AEDC dump combustor simulation facility. 
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ERROR RESIDUALS 
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Fig. 16 Convergence history for the non-reacting 
dump combustor calculation. 
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Fig. 17 Comparison of computed and experimental 
centerline velocities. □ , experiment; solid line, 
k-«{Tf model; dashed line, standard k-« model. 
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the centerline velocities well in the initial region (X/D=0 to 
X/D=l) and the final region (X/D=5 to X/D=6). Near the measured 
reattachment point (X/D=3), the standard k-e model underpredicted 
the velocity and the k-e{T) model overpredicted the velocity. 

The measured reattachment length for the non-reactive dump 
combustor corresponds to a step height of 9.9 at an expansion 
ratio of 2:5. While Chaturvedi's (ref. 37) pipe expansion 
experiment has measured a 9.2 step heights for an expansion ratio 
of 1:2. The standard k-e turbulence predicted a shorter 
reattachment length at 8.32 step heights, and the k-e(T) 
turbulence model predicted a longer reattachment length at 10.7 
step heights in this study. The longer measured reattachment 
length for the non-reactive dump combustor is certainly 
reasonable for a larger expansion ratio, in comparison to 
Chaturvedi's measurement. In general, predictions made by both 
models were within the error bounds of the experimental data. 

Fig. 18 gives the comparison of the axial distributions of 
wall static pressures. The predictions agreed well with the 
experimental data at both ends. From X/D®1 to X/D=4, the 
standard k-e model overpredicted the pressure and the k-e(T) 
model underpredicted the pressure. 

Fig. 19 gives the comparison of the axial distribution of 
hydrogen mass fraction along the wall. Both turbulence models 
predicted complete mixing after four diameters downstream. The 
experiment and prediction agreed reasonably well beyond one 
diameter downstream. The suddenly drop of the hydrogen mass 
fraction near about 0.5 diameter downstream is believed to have 
caused by the formation of the recirculation zone. 

A comparison of the radial distribution of mean axial 
velocity at several axial stations downstream of the step is 
given in Fig. 20. Both turbulence models predicted good results 
in comparison with the experimental data. The predicted 
reattachment point was about three diameters downstream of the 
step that matched the experimental data. Comparison of the 
computed and experimental radial profile of hydrogen mass 
fraction are shown in Fig. 21. The predictions agreed well with 
the data. Although there was a small difference of about 0.001 
between the predicted values and the experimental values at the 
last station. The predicted values are well within the error 
bound of the experiment. 
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Fig. 18 Axial distribution of wall static pressure, 
o, experiment; solid line. k-«|T| model; 
dashed line, standard k-« model. 
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Fif. 19 Axial decay of hydrogen mass fraction at wall. 
O. experiment; solid line. k-«}T| model; 
dashed line, standard k-« model. 
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Fig. 20 Radial distribution of mean axial velocity. 
□ , experiment; solid line, k-«|T| model; 
dashed line, standard k-« model. 
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Fig. 21 Radial distribution of hydrogen mass fraction. 
O, experiment; solid line, k-«|T| model; 
dashed line, standard k-« model. 
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A Reactive Hvdroqen/Air Ramjet Dump Combustor 


Results of the Reactive Dump Combustor Calculation 

The experimental results of the reactive dump combustor have 
shown a longer recirculation zone due to the heat release of the 
exothermic chemical reaction. The computational domain of the 
reactive dump combustor calculation was therefore extended to 
eight diameters downstream of the step so that the longer 
recirculation zone could be properly contained. The typical grid 
size used was 70x31. A grid size of 105x31 to include fifteen 
diameters downstream of the step was also computed, but the 
results showed little or no difference with the 70x31 case. 

Fig. 22 shows the comparison of computed and experimental 
centerline velocities for the reactive dump combustor 
calculation. The k-e{T) prediction gives excellent comparison 
with the data. It was used exclusively for latter comparisons. 
The predictions with the standard k-e and the extended k-e 
turbulence models showed too rapid a mixing rate without 
considering the temperature effect. The effect of the k-e{T} 
turbulence was to reduce the turbulent eddy viscosity. Fig. 23 
gives the comparison of axial wall static pressures. The 
prediction and the experimental data are in good agreement. 

The comparison of experimental and computed radial 
distributions of mean axial velocity are shown in Fig. 24. The 
computed values are in good agreement with the experimental data. 
The reattachment length predicted was a little less than six 
diameters after the step and matched the reported data. In 
comparison, Drummond's predictions showed no reversed axial 
velocity at any the axial stations. The k-e{T) turbulence is 
obviously a better representation for the chemically reacting, 
recirculating flow calculation. Fig. 25 gives the comparison of 
the computed and experimental radial temperature profiles. 
Excellent agreement was obtained between the predicted values and 
the experimental data at all stations. The predicted flame front 
was located approximately between X/D=l and X/D=2. The heat 
release of the reaction was therefore properly represented by the 
one-step reversible chemical kinetics model. In comparison, 
Drummond's predictions show an early radial temperature rise and 
then a drop of temperature as it approaches the wall at X/D-l, 
and the wall temperatures were consistently lower than the 
experimental value for the rest of the stations. This was 
probably caused by a combined effect of the chemistry and 
turbulence models. The complete combustion chemistry forced a 
flame front very close to the step, and the algebraic turbulence 
model produced fast mixing that did not allow enough unreacted 
fuel to remain near the wall. The predicted hydrogen mass 
fraction profiles were in reasonable agreement with measured 
data. The required CPU time for the reactive dump combustor 
calculation was approximately three times longer than that of the 
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non-reactive dump combustor calculation. The convergence 
criterion was relaxed to three orders of magnitude decrease in 
the residuals since further integration did not improve the 
results. 

In comparison of the non-reactive and the reactive cases, 
both the experiments and the predictions indicated that the heat 
release of the chemical reaction have reduced the effective eddy 
viscosity and increased the length of the recirculation length 
for the dump combustor configurations. 
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Fi£. 22 Comparison of computed and experimental 
centerline velocities, o, experiment; solid line, 
k-«|T| model; dashed line, extended k-« model; 
dotted line, standard k-« model. 
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Fig. 23 Axial distribution of wall static pressure. 
O, experiment; solid line, k— « |T| model. 
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Fig. 24 Radial distribution of mean axial 
velocity for reactive dump combustor. □ , 
experiment; solid line, k-«|Tj model; dashed 
line, Drummond’s calculation. 
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Fig. 25 Radial distribution of mean static 
temperature. □ , experiment; solid line, 
k-«{Tj model; dashed line, Drummond's 
calculation. 
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TRANSIENT FUEL PREBURNER SIMULATION 


Transient Upstream Boundary Conditions 

At the start command, the initial flowfield of the fuel 
preburner is filled with atmospheric nitrogen purge gas as 
described by ref. 38. After the start command, the fuel and 
oxidizer flow out of the augmented spark igniter (ASI) and the 
injector elements on the faceplate, to prime the chamber and the 
turbine. The transient upstream boundary conditions therefore 
include the ASI and injector fuel and oxidizer mass flow rates, 
and the chamber pressure. These transient flow properties were 
obtained from Digital Transient Model (ref. 39) simulation (ref. 
8). They are written as: 


P 

W 

W 


ASI, i 
INJ, i 


P r (t) 


W 


ASI , i 


(t) 


W INJ, i < 


(31) 


where subscript i represents 
The species mass fractions at the 


either fuel(f) or oxidizer(o) . 
inlet boundary are given by 


«i = Wi/IWj (32) 

and the total mass flow rate at the inlet is 


W T = IW* (33) 

The operational pressure variable p in FDNS2DR is the 
difference between the actual pressure and the pressure at a 
reference point. When the reference point is placed at the inlet 
boundary and the reference pressure is varied with the transient 
chamber pressure, the actual pressure in the whole flowfield will 
respond instantaneously to the transient chamber pressure and to 
the solution of the conservation equations. The density at the 
inlet boundary is then calculated as: 


p' = (P + PcJ/tZai/Mi RT) (34) 

while the velocities at the inlet are updated by 


u' = uW T //p'udS 

and (35) 


v ' = 0 


Fig. 26-28 show the DTM simulated transient flow properties 
for a nominal operation at an elapsed time of one second. The 
engine main fuel valve begins to ramp open at the start command, 
and the FPB oxidizer valve begins to open at 0.1 seconds (ref. 
3). The ASI flow rates are shown in Fig. 26. The fuel flow did 
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Fig. 26 Transient fuel preburner ASI inlet mass flow rate 
boundary conditions. 
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Fig. 27 Transient fuel preburner chamber pressure boundary 
condition. 
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Fig. 28 Transient fuel preburner injector mass flow rate 
boundary conditions. 
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not emerge until 0.07 seconds and the oxidizer flow did not 
appear until 0.15 seconds. The oxidizer flow rate had a drop 
from 0.4 to 0.6 seconds, and the fuel flow rate oscillated at 
this time period. The chamber pressure history is given in Fig. 
27. The pressure did not start to rise until 0.45 seconds. Fig. 
28 shows the DTM simulated injector flow rates. There is an 
initial decline of fuel flow rate at around 0.45 seconds, 
corresponding to the time of pressure rise in Fig. 27. A major 
second drop of fuel flow rate occurred near 0.6 to 0.8 seconds. 

In the meantime, an increase of oxidizer flow rate also occurred. 
It can be seen that the fuel flow rates were probably designed to 
have a lead over the oxidizer flow rates to avoid a local lean 
mixture in the chamber. But the pressure wave sent from 
downstream of the fuel preburner due to the sudden pressure rise 
may have shut down the fuel supply line for a short period of 
time. During this period of time, the fuel flow rate dropped, 
and, the oxidizer flow rate increased in order to compensate for 
the fuel flow rate loss. The combination of the fuel flow rate 
drop and the oxidizer flow rate increase have caused a lean 
mixture. The lean mixture combustion would cause a temperature 
spike. The closer the mixture ratio is to the stoichiometric 
composition, the higher the thermal gradient will be. 


Grid Definition and flow assumptions 

The computational grid was constructed based on the 
assumption that two concentric jets can be used to represent the 
fuel preburner flows. The inner jet represents the ASI flow and 
the outer jet stands for the injector flow. Gaseous propellants 
were assumed to come out of the injectors on the faceplate. Many 
grids were designed and evaluated with the technique of 
quantitative measures (ref. 40) and series of steady state flow 
tests, so that the final computational grid can withstand the 
transient, swirl, reactive flow with temperature spikes. Fig. 29 
shows the final computational grid for the calculation. The grid 
size was 68x45 and no grid related problems has occurred during 
the calculation. 

FDNS2DR methodology was used to describe the response of an 
axisymmetric FPB and turbine inlet thermal flowfield to the 
transient start. To simulate start-up, the propellants must be 
ignited by the ASI system and the flame must propagate through 
the chamber at a realistic rate as determined by the chemical 
kinetics and mixedness level of the system. The upstream 
boundary conditions were obtained from the DTM simulation, and 
therefore simulate the actual opening of the FPB0V and the 
response of the transient combusting flow. Downstream boundary 
conditions at the FPB exit are extrapolated from the interior 
flowfield and adjusted by the transient inlet mass flow rate. 

This does represent an additional degree of approximation but it 
was found necessary to obtain a well behaved solution. Adiabatic 
walls were also assumed for the entire simulation. 
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Results of the Transient Calculation 


The transient fuel preburner thermal flowfield calculation 
was carried out for 1.0 second. Fig. 30 - 40 show the calculated 
FPB thermal flowfield time histories at an interval of 0.1 
seconds. Fig. 30 shows the thermal flowfield at the start 
command. The upper half of the preburner represents the flow 
vector streak lines, while the lower half gives the temperature 
contours. The left-hand side of the preburner indicates the 
inlet of the propellants. The circular area represents the top 
of the turbopump dome (Kaiser hat). It can be seen from the 
temperature contour that the hot igniter flow was coming out of 
the AS I chamber, and the streak lines showed a quiescent nitrogen 
gas filled environment. Fig. 31 shows that the igniter flow 
flame was swirling out of the ASI chamber and penetrated into the 
injector flow region, at about 0.1 seconds. The flow in the rest 
of the chamber was still undisturbed nitrogen purge gas. At 0.2 
seconds. Fig. 32 shows the hot igniter flow was advancing farther 
downstream. The injector flow appeared to have gained some 
momentum, but was not combusted since it was mostly cold 
hydrogen. About half of the nitrogen purge gas was pushed out of 
the combustor at this moment. At 0.3 seconds, Fig. 33 shows the 
hot ASI flame has extended to the top of the turbopump dome. 

Most of the nitrogen gas was pushed out. The injector flow 
gained more momentum but was still cold since the injector 
oxidizer had just entered the chamber. Fig. 34 shows the igniter 
flow had covered the whole Kaiser hat and a flame front had begun 
to move toward the faceplate. The injector oxidizer flow rate 
had increased. The nitrogen concentration had dropped to a 
negligible amount and the species nitrogen was turned off for 
more efficient calculation at 0.4 seconds. 

Fig. 35 shows the thermal flowfield corresponding to about 
0.5 seconds after the start command. The initial decline of 
injector fuel supply had made the mixture ratio leaner. The 
flame front propagated to the faceplate. A hotter region had 
formed on top of the turbopump dome, and propagated toward the 
ASI. The flow had sped up due to the increase in flow rate and 
the combustion. Notice a centerline recirculation zone was being 
formed near the nut location of the Kaiser hat. The calculation 
proceeded smoothly although the DTM inlet igniter fuel flow rate 
had shown some oscillations due to the sudden increase of the 
chamber pressure. At 0.6 seconds, the injector fuel flow rate 
had past the initial decline and was restored to its original 
level. Notice the fuel rich combustion and the sudden increase 
in chamber pressure had slowed down the flow vectors in the 
chamber. The hotter spot of the igniter flow had continued to 
propagate towards the ASI, as shown in Fig. 36. 

At 0.6 seconds after time zero, the fuel flow rate begins 
its second drop to a mixture equivalence ratio of about five. It 
would rise a little bit but proceed to drop back to the same 
equivalence ratio of about five at 0.7 seconds. Fig. 37 shows 
the thermal flowfield at 0.7 seconds. The igniter flow flame 
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Fig. 29 Finite difference grid used to compute transient fuel 
preburner flowfield. 
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Fig. 30. 



Fig. 31. 
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Fig. 32. 



Fig. 33. 
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Fig. 37. 
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front had extended all the way back to the exit of the ASI 
chamber. Hot spots formed near the Kaiser hat nut, and hot 
streaks flowed over the turbopump dome. The hottest temperature 
(2700 deg. R) matches the DTM calculated peak temperature at 0.7 
seconds for the nominal start. A disastrous start would have 
Introduced near stoichiometric combustion in the combustion 
chamber. After 0.8 seconds, fuel rich combustion restored as a 
consequence of the Increase in fuel flow rate. This is shown in 
Fig. 38.-40. In the present axisymmetric FPB study, the nut of 
the Kaiser hat seemed to withstand a prolonged period of high 
thermal gradients . 

The transient radial temperature distribution at the fuel 
preburner exit is given in Fig. 41. The temperature rose 500 
deg. R at about 0.5 seconds. A maximum temperature is reached at 
about 0.7 seconds. The radial temperature dropped afterwards due 
to the increase in fuel flow rate. Notice the temperature is 
hotter near the Kaiser hat wall side. This implies that the 
turbine blade fillets will undergo higher thermal gradients. 

Fig. 42 represents the transient temperature profile along the 
Kaiser hat wall. The hottest temperature occurred near the 
Kaiser hat nut at about 0.8 seconds. The centerline temperature 
distribution is shown in Fig. 43. A flame front started near the 
Kaiser-hat nut and moved toward the igniter flow inlet as time 
went on. The flame was then quenched after 0.8 seconds. The 
overall centerline temperature was near its hottest at about 0.7 
seconds. The corresponding ASI exit temperature had reached 
about 1700 deg. R. The maximum temperature occurred near the nut 
at 0.8 seconds. This is probably due to the centerline 
recirculation zone formed in front of the nut. 

A comparison of the FPB start-up transient temperature 
history at the fuel preburner exit and turbine inlet is shown in 
Fig. 44. Calculated temperature transients are given at two 
locations: at the middle of the exit line, and at the exit on the 
Kaiser-hat wall. The measured (ref. 3) turbine inlet temperature 
is represented by the thermocouple located near the inboard side 
of the turbine inlet and in between the baffles, so that three 
dimensional effects are minimized. Both the calculation and the 
measurement showed a peak at about 0 . 5 seconds , and a second peak 
at about 0.7 seconds. The temperature was then greatly reduced 
as the time approached one second. The measured first peak 
temperature is lower than that of the calculations. The final 
peak temperature and the peaking time of the FDNS2DR calculation 
at the middle of the FPB exit agreed very well with those of the 
measurement. The exit temperature calculated on the Kaiser hat 
wall is higher than that calculated at the middle of the exit, 
since the hotter striated torch flow was flowing near the 
turbopump dome. In general, the transient temperature history of 
the calculation showed qualitatively excellent agreement with 
that of the measurement . 

The transient fuel preburner species mass fraction history 
at Kaiser hat nut is given in Fig. 45. The unburned oxygen mass 
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fraction dropped at 0.5 seconds, indicating the onset of the 
combustion. The mass fraction curve of the combustion product, 
H20, peaked at about 0.78 seconds, representing the highest 
temperature in the time frame considered. The unburned hydrogen 
mass fraction also dropped to its local minimum at 0.78 seconds. 
These three species mass fraction curves indicate that the 
combustion reaction near the Kaiser hat nut was almost complete, 
possibly with the help of flow recirculation. Similar curves are 
shown in Fig. 46 for the species mass fractions at the center of 
the fuel preburner exit. H20 mass fraction peaked and unburned 
hydrogen bottomed at an earlier 0.73 seconds. The unburned 
oxygen curve indicated an incomplete combustion, since ten to 
fifteen percent of oxygen was left unburned after 0.45 seconds. 
This is probably caused by the overall lower temperature in the 
injector region. Since the calculated temperature history at the 
fuel preburner exit based on the species composition was very 
closed to the measured turbine inlet temperature history, as 
previously shown in Fig. 44. The calculated species composition 
history is believed to be very close to the actual transient 
operation. Fig. 45 and 46 showed two different types of 
transient species composition history at two different flow 
situations. Obviously the remaining question is whether the 
unburned oxygen at the exit will cause near stoichiometric 
combustion in the complex flowfield of turbine stages, if so, 
this is certainly a likely source of initial turbine blade 
cracking and should be carefully studied with reactive CFD 
method . 
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Fig. 41 Simulated radial temperature history at the fuel 
preburner exit. 
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Fif. 43 Simulated centerline temperature history from ASI 
exit to Kaiser hat nut. 
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Fig. 44 Comparison of temperature history at calculated fuel 
preburner exit and measured turbine inlet. 
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Fig. 45 Transient fuel preburner species mass fraction 
history at Kaiser hat nut. 
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Fig. 46 Transient fuel preburner species mass fraction 
bistory at middle of exit. 
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Conclusions and Recommendations 


A pressure based finite rate, reactive CFD model FDNS2DR has 
been developed to simulate the start transient of the SSME FPB 
operation. An axisymmetric configuration with a grid size of 
68x46 was used to represent the FPB flowfield. The transient 
upstream boundary conditions were obtained from a one-dimensional 
DTM simulation to better represent the actual operation. The 
results of the CFD calculation have shown temperature spikes near 
the FPB exit. The timing and magnitude of the temperature spikes 
agreed well with those of the Aerojet measured turbine inlet 
temperature data. In addition, an appreciable amount of oxidizer 
was left unburned to enter the turbine stages. This results can 
not be obtained from a steady state CFD thermal flowfield 
simulation. Therefore, FDNS2DR is a useful CFD engineering 
design tool for studying the transient operation dominated SSME 
flow anomalies. 

As a result of the transient thermal flowfields described by 
the analyses reported herein, the following future work is 
recommended : 

1. Simulate a full three dimensional transient fuel preburner 
calculation so that the effect of baffle hydrogen flow can be 
included. 

2. Parametrically run the DTM simulation so that the effect of 
different oxidizer valve sequences can be studied in conjunction 
with the downstream pressure wave characteristics. Various 
numerical experiments can also be performed with FDNS2DR to alter 
the ignition and shut down characteristics to study possible 
modifications to the operating procedure in order to attenuate 
the thermal shocks. 

3. Although the thermal gradients in shutdown transient are less 
severe, but possibly more significant since the duration is much 
longer. A shutdown transient study is recommended for the 
complete under standing of the turbine blade crack problems. 

4. Study the combustion characteristics of the unburned oxidizer 
once it enters the turbine stages. 
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THERMODYNAMIC AND TRANSPORT PROPERTIES FOR OXYGEN 


• Introduction 

Thermodynamic properties of single-component fluids are 
defined in terms of thermal and caloric equations of state. At low 
pressures or high molal volumes, all substances are gases which 
obey the perfect eras equation of state. 

P = RT/MV 

where P is pressure, R is the universal gas constant, T is 
temperature, M is molecular mass, and V is specific volume 
(1/density) . The caloric equation of state relates enthalpy, H, 
or internal energy, U, to two state variables, say temperature and 
density. If the internal energy is a function of temperature only, 
the gas is said to be ideal . If the gas is thermally perfect and 
calorically ideal, 

Cp/C v = 1 + R/C v 

where C p and C v are constant pressure and constant volume heat 
capacities, respectively. If the heat capacities are constant, the 
ratio 


7 = Cp/C v 

is conveniently used to represent all of the thermodynamic 
properties of the gas. In general, as molal volumes decrease or 
density increase these simple relationships become invalid; at even 
higher pressures the "gases" condense to form liquids. The theorem 
of corresponding states postulates that the critical conditions may 
be used to non-dimensionalize the equations of state so that all 
pure gases exhibit the same compressibility factors at the same 
reduced density and temperature. The compressibility is a 
multiplitive term included in the perfect gas law to make it apply 
to any temperature and density. Figures 1 and 2 show a non- 
dimensionalized plot of 0 2 properties as predicted by the methods 
described herein which are in agreement with the data reported by 
NBS (1,2). Using a third parameter, the critical compressibility, 
pressure and enthalpy are well correlated in terms of reduced 
density and temperature for many substances (3) . A more convenient 
correlation is that reported in (4,5) as the HBMS equation of 
state, since all of the correlations are expressed as analytic 
functions. Representations of the NBS thermodynamic data with the 
HBMS model and of transport property empirical correlation 
equations will be summarized so that the required properties for 
two-phase flow analysis of oxygen can be provided for a numerical 
solution to the conservation equations. 
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• The Thermal Equation of State 

The pressure data shown in Figure 1 may be divided into a gas 
region (I) , a dense gas region (II) , a liquid region (III) , and a 
two-phase region (IV) . Solid oxygen will not be considered. The 
gas region may be represented with a modification to the perfect 
gas law. Van der Waals ' equation of state accounts for the energy 
of attraction between molecules with the "a" term and for physical 
volume occupied by the molecules with the "b" term. 

(P + a/V 2 ) (V-b) = RT 

The HBMS equation generalizes this equation, thusly 

[P + a{T}/V 2 + a' (T)/V 3 ) [V - b + b'/V] « rt 

Temperature, density, and pressure are non-dimensionalized by 
dividing by the critical temperature, density and pressure, 
respectively. These variables are denoted by t, p, and p, 
respectively. Braces are used herein to denote functionality. The 
functions a and a 1 , and the constants b and b' were evaluated in 
(4) to produce the HBMS/I equation of state. 

Pi = -(ko+(/3-ko)t" l )p 2 - O.5(l-k o -ct+20) (t+t _1 )p 3 + 

(1+/?) 3 pt/ (0(30-1) - (30 3 -60-l) p+0 (0-3) p 2 ] 

where k 0 is 5.5, a = 7 +(z e _1 - 3.72)/0.26, and 0 is defined by 

Z c = P c V c /RT e = 0(30-1)/ (1+0) 3 

Alternative expressions are given for k 0 and a in (4) , but they 
are not needed to represent 0 2 properties. 

The expression for pressure required modification to represent 
the dense gas region, so HBMS/II was developed (4) . 

Pn ~ t 1 [B 30 p 2 + B il0 p 3 ] + 

[BdP* 1 + B u + B 21 p + B 31 p 2 + B 41 p 3 ] + 

t[B„ 2 p _1 + B 12 + B 22 p + B 32 p 2 + B A2 p 3 + B 52 p*] + 

t 2 [B 03 p _1 + B 13 + B 23 p + B 33 p 2 ] 

In HBMS's nomenclature, this equation is 

Pn = E t j_1 Qj { p } where j « 0, 1, 2, 3. 

For the liquid region, HBMS/III is 

Pm = 2 t J_1 (Qj{p ) - Qj{p l >) + Pv where j - 0, 1, 2, 3. 
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The saturated liquid density, p L , and the vapor pressure, Pv, must 
be obtained from additional correlation equations. Riedel's vapor 
pressure equation was used in (4) because of its accuracy near the 
critical point. 

In {pv} = a In {t} + 0.0838(a - 3.75) (set* 1 - 35 - t 6 + 

42 In { t ) ) 

This equation becomes progressively less accurate as the normal 
boiling point is approached. Attempts to interpret a as a 
parameter to fit the vapor pressure at arbitrary temperatures 
destroy the accuracy of the fit near the critical point. other 
vapor pressure correlations should be considered if the region of 
interest includes pressures which are less than 0.4 T c ? the NBS 
report (1) suggests a high order polynominal which is used herein. 

In {P} = K x + A 2 T + A 3 T 2 + A 4 T 3 + A s T 4 + A 6 T 5 + A 7 T 6 + A 8 T 7 
where P is in atm and T in Kelvins. The A's are: 


A x = -62.5967185 
A 3 = -4 . 68973315E-2 
A s = -4 . 09349868E-6 
A 7 = -5 . 13113688E-11 


A 2 = 2.47450429 
A 4 = 5 . 48202337E-4 
A 6 = 1 . 91471914E-8 
A e = 6. 02656934E-14 


Predictions with these two vapor pressure correlations are 
shown in Figure 3. The polynominal equation fits the data over a 
wider range of temperatures. If an accurate vapor pressure 
correlation is not available for a specific gas of interest, such 
a correlation should be developed. 


The saturated liquid density may be represented by 
P L {t) = 1 + c(l-t) 1/3 + d(l-t) 


where c and d are evaluated by using two data points from the 
saturated liquid curve (1) . Predictions with this equation are 
compared to the NBS data in Figure 4; the comparison is very good. 

The working form of the HBMS/III pressure equation is 

Piii “ t 1 [B3oP 2 + B 40 p 3 ] + 

[B 01 p 1 + B u + B 21 p + B 31 p 2 + B 41 p 3 ] + 

t[B 02 p 1 + B 12 + B 22 p + B 32 p 2 + B 42 p 3 + B 52 p 4 ] + 

t 2 tB 03 p' 1 + B 13 + B 23 p + B„p 2 ] - 

t' l [B 30 p L 2 + B 40 p L 3 ] - 
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[ B oiPl 1 + ®11 + ®2lPL + ®3lPL 2 + B uPl 3 ] “ 

^[602^1. 1 + B 12 + B 22PL + B 32 PL 2 + B «PL 3 + B 52 Pl *3 “ 

^ [ B 03Pl B 13 ■*" ®23Pl B 33Pl 2 ] ^ Pv(t') 

The values of the B^'s are given below. These values 
correspond to those given in the text of (4) ? they do not match the 
values given in the tables of either (4) or (5) which apparently 
contain some algebraic errors. 


B 30 

— 

5.5 - p 


B *0 

s= 

-2.25 - 0.5 a + 0 


B 01 

= 

88.5 - 3.12 p 


B n 

= 

-313.3 + 13.42 P 


B 21 

= 

408.9 - 21.54 P 


B 31 

= 

-237.4 + 15.30 /3 


B« 

= 

47.8 - 4.06 p 


B 02 

= 

-124.46 + 3.84 p + 0.363 p 2 

B 12 

= 

429.0 - 9.84 P - 1.815 

P 2 

B 22 

= 

-528.8 + 1.98 P + 0.363 p 2 

B 32 

= 

263.0 + 15.70 p - 3.63 

P 2 

B 42 

= 

-27.05 + 0.5 a - 16.18 

P + 1.815 p 

B 52 

= 

-8.44 + 4.50 P - 0.363 

P 2 

B 0 3 

= 

44.4 - 5.22 p 


B 13 

= 

-156.9 + 18.92 p 


B 2 3 

= 

204.3 - 25.44 p 


B 33 

= 

-115.5 + 15.0 p 


B 43 

= 

23.7 - 3.2 p 



Predicted pressures as a function of temperature and density 
are shown in Figure 1. The HBMS thermal eguation of state is an 
excellent representation of the NBS data. 
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• The Caloric Equation of State 

For the gas region, the internal energy excess function may 
be written as 

(U - U 0 )/RT = z e J' [ £ - ( M j j % 

Since H = U + PV, 

(H x - H 0 )/RT = Zc J P [ t ” ( ft ) ] J* + (Z ° P/Pt) ” 1 

For the dense gas region, the enthalpy excess function is 

(Hu - H 0 )/RT = z c J' J j j % + ic jE " Pi j 

+ (H x - H 0 )/RT 

where p x and H x are evaluated at the intersection between the gas 
and the dense gas regions. For the liquid region, 

(H m - H 0 )/RT * Z C J P [t”[ft] ] ^ + ^[ E_Ev ] 

+ (H l - H 0 )/RT 

Inserting the appropriate pressure equation in the enthalpy 
excess function, the following working equations are obtained. 

(H x - H 0 )/RT = - Z e [2k 0 t' 1 + 3 (/3-k 0 ) t‘ 2 ) p + 0. 5 (l-k 0 -a+2/S) (l-2t' z ) p 2 ] 

+ (l-bp+b'p 2 )" 1 + 1 

where b = (3/3 2 -6/3-l)//3 (30-1) and b* = (/3-3) (3/3-1) . 

(H n - H 0 )/RT = t~ 2 [ 3B 30 p + 2B, 0 p 2 ) + t'^O. 5B 0lP ‘ 2 + B 21 In (p) + 2B 3X , 

+ 1.5B 41 p 2 ] + [B 02 p' 2 + B 12 p _1 + B 32 + B„p 2 + B 52 p 3 ] - 
t' 2 [ 3B 30 + 2B 40 ] - f^O-SBo: + 2B 31 + 1.5B 41 ) - [B 02 + 
b 12 + b 32 + b* 2 + B 32 ] + Z e t‘ 1 tp II p" 1 - Pi] + 

(H 3 - H 0 )/RT 

where p 3 = -(ko+^-koJt* 1 ) - O.5(l-ko-or+20) (t+t' 1 ) + 

( 1+0) 3 t/ [0 ( 3/3-1) - (30 3 -60-l) +fi (/3-3 ) ] 

and 
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(Hi - H 0 )/RT = - z c [ (2k 0 t' 1 + 3(^-ko)t' 2 ) + O.5(l-)c 0 -a+2j8) (l-2t' z ) ] 

+ (l-b+b')' 1 + 1 
For the liquid region, 

(H m - H 0 ) /RT = 2 c [t' 2 [3B 3oP + 2B< 0 p 2 ) + t _1 [ 0 . 5B 01 p' 2 + B 21 In <p} + 2B 31 , 
+ 1.5B 41 p 2 ] + [B 02 p* z + B 12 p _1 + B 32 + B* 2 p 2 + B 52 p 3 ] 

- t' 2 [3B 30 p L + 2B, 0 p L 2 ] - t" 1 [ 0 . 5B 01 p L ” 2 + B 21 In {p L } + 


2B 31 p L + 1 • 5B 41 p L *] — [B 02 p L 2 + B 12 p L + B 32 + B„p L + 


2i 

L J 

-”2 i 


B 52 p l 3 ]] - Z c [2f 2 (B 30 p L 2 + B A0 p L 3 ) + t* 1 (B 01 p L _1 + B n + 
®2iPL + B 31 p L 2 + B u p L 3 ) - t{B 03 p L * 1 + B 13 + B 23 p L + 


B 33 P L Z + B* 3 p L B )](p t " 1 - p' 1 ) + Z e [t’*(2B 30 p L + 3B w p t *) + 


-l 


(“B 0 iPl" 2 + B 21 + 2B 31 p L + 3B 41 p L 3 ) + t(-B 02 Pl “ + B 22 + 
2B 32 Pl + 3B* 2 p L + 4B 52 p L ) + t ( w B 03 p L + B 23 + 2B 33 p L + 
3B 43 p L 2 ) ] (p L _1 - P* 1 ) JsUlJ + Z e (p m - Pv)/pt 

- Z c ( P l " 1 - p" 1 ) gEvj 

This equation may be algebraically simplified, but the form 
presented expedites eliminating computational errors. To use this 
equation, the temperature derivatives of p L and py and H L are 
required. The p L derivative is accurately described by the 
previously stated correlation equation. H L may be calculated from 
the Clapeyron equation 

(H v - H l )/RT = z c (p/ 1 - p^JgBvj 

Riedel's vapor pressure equation is not sufficiently accurate to 
calculate the indicated derivative at low temperatures. The 
polynominal equation is sufficiently accurate for the enthalpy 
correction, but not for evaluating H L near the critical point. The 
scheme recommended is to use Riedel's equation for evaluating H L at 
a t of 0.9 and to then use the correlation from (6) to determine 
H l at any other temperature. Specifically 


-2 


(H v " - H l ) x 


(1 , O-T - t zl 
[l.O - tj 


0.38 


where the subscript 1 indicates the value at t = 0.9 and the 
subscript 2 indicates the value at the desired temperature. 

The predicted heat of vaporization is compared to the NBS values 
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in Figure 5. This correlation for H* may be used with the 
Clopeyron equation to develop a correlation for P v between t = 0.9 
and the normal boiling point. 

The final enthalpy predictions are compare to NBS enthalpy 
data in Figure 2? this comparison is very good. The dimensionless 
excess function is compared to tabulated values of this parameter 
from (3) in Figure 6. The shape of the two correlations is very 
close, but some differences are seen in the magnitude of the 
correction. The comparison in Figure 2 indicates that the HBMS 
equation of state yields better enthalpy values than the 
correlations from (3) . 


• The Two-Phase Region 

The quality, s, and saturation properties are related by 
P iv — [s/pv + (1.0 - s )/p L ] 1 
H IV = H v s + H l ( 1.0 - s) 

If p and t or H are specified, the state is determined. Specifying 
p and t (or H) does not uniquely determine the state in the two- 
phase region. 


• Low Pressure Enthalpy 

H 0 is the enthalpy at low pressure which is that of an ideal 
gas, hence 



Many heat capacity expressions have been reported which use a lower 
limit of a non-zero reference temperature (including the NBS data 
(1)); such expressions are not convenient for general flow 
calculations. The heat capacity correlation recommended is that 
from (6) 


C p = C pa + C pb T + C pc T 2 + Cp,, T 3 
The units on C p are Btu/lb n R, and the constants are 


C pa = 0.20978 
C pb = -7 . 6040E-9 
C pe = 1 . 3407E-8 
Cpd = -3 . 4079E-12 
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• Transport Properties 

The thermal conductivity is represented by three additive 

terms 


k = ko{T} + k.{p,T} + k c { p , T } 

where k„ is the low pressure gas contribution, k. is the excess or 
dense gas contribution, and k c is due to enhancement near the 
critical point. Values for the thermal conductivity taken from (1) 
are shown in Figure 7. Analytical expressions for k„ and k, are 
given in reference (7) . These expressions are: 

ko = J Xj T 1 for i * 1, * . . , 9 

k, s S Xj P i for i = 10,.., 15 

A convenient expression was not given for k c , therefore step 
function about the critical point is recommended for describing 
superconductivity. An 8-fold increase was assumed for the critical 
point; this factor decreases to zero 2 degrees away from the 
critical temperature. This increase was estimated from C0 2 data 
given in (6) . The correction is used only at these temperatures 
and if the density is within ± 5% of the critical value. 

The viscosity does not experience abnormal behavior near the 
critical point so only two terms are needed, one each for the low 
pressure gas and for the dense gas. 

M = Mo(T} + M.(p, T} 

Analytical expressions are given for these terms in (7) ; values for 
the viscosity are given in Figure 8. Scant data exist for bulk 
viscosity, therefore assume it is zero and use only the viscosity 
to evaluate normal stresses. Normal stresses will be important 
only for low density gas flows where it is known that the bulk 
viscosity is small. Note that an erroneous definition is used for 
bulk viscosity in the NBS report (1), which implies that assuming 
zero bulk viscosity eliminates all normal stresses. Such is not 
the case. The correlations for viscosity are: 

M 0 = £ Yi T 1 for i = 1 , . . . , 9 

M. = 2 y A p i for i = 10,.., 12 

Notice that temperature effects on the dense gas contributions 
to thermal conductivity and viscosity do not appear in the 
correlations, such effects have not yet been established. The 
constants are: 

X x = 9.3648096074E-7 y x = 8 . 344128902E-7 

X 2 = -4.1398868586E-9 y 2 = -6. 4319584704E-9 
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1 . 338584 04 12E-10 


*3 = 

1 . 02783713 25E-10 

X, = 

-1 . 059 07 3 04 3 4E- 12 

*5 = 

5 . 9317072836E-15 

*6 = 

-1.9858382126E-17 

II 

r> 

X 

3 . 98094 22 422E-20 

X 8 = 

-4.4173352317E-23 

Xg = 

2 . 0893500738E-26 

II 

O 

H 

X 

6 . 2807950678E-4 

*11 = 

-4 . 9337311644E-4 

X :2 = 

2 . 5242925212E-3 

X 

»-» 

o> 

II 

-5. 1527923868E-3 

X 14 = 

5 . 4461096782E-3 

Xl5 = 

—1 . 8991126599E-3 


y 3 = 

y* = -1.3283970709E-12 
y 5 - 7.4604203289E-15 
y 6 * -2.5398419748E-17 
y 7 = 5.2077731498E-20 
y B = -5.9258985472E-23 
y 9 = 2.8773126492E-26 
y 10 = 4 . 7293376329E-4 
y :i « -1 . 7410413651E-4 
y 12 - 5.9995361171E-4 
y 13 = 0.65391907848 
y u = 2.9886313449E-5 


y 15 = 9.25 

The units on k are W/cm K, and on n are gm/cm s. T is in K, 
and p is in gm/cm 3 . If the density is above 0.932 gm/cm 3 , a 
different correlation equation is used: 


Me = (y i3 p + yu(exp{y 15 ,}-1.0))/1000 

If different units are desired, it is more convenient to make the 
change once the transport property has been calculated. The 
correlation equations give transport property values comparable to 
those shown in Figures 7 and 8. 

Transport properties in the two-phase region are determined 
from the quality with an equation like that used to determine 
enthalpy. 
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♦ Closure 


The property equations presented herein have been evaluated 
with a FORTRAN program. The formulation has been presented so that 
a minimum of terms have to be recomputed during the course of 
evaluating all of the properties and so that the empirical 
constants can be identified in terms of the original work from 
which they were obtained. The properties estimated with the 
aforementioned code are not intended to be more accurate than the 
referenced NBS values, but they are of good engineering accuracy 
and are of a form which is convenient to use when computing 
flowfields. The FORTRAN code produces files of properties which 
may be used as tables of data for a CFD calculation. 
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